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Abstract 

A classical problem in acoustic scattering concerns the evaluation of the Green's function for the 
Helmholtz equation subject to impedance boundary conditions on a half-space. The two principal ap- 
proaches used for representing this Green's function are the Sommerfeld integral and the (closely related) 
method of complex images. The former is extremely efficient when the source is at some distance from 
the half-space boundary, but involves an unwieldy range of integration as the source gets closer and 
closer. Complex image-based methods, on the other hand, can be quite efficient when the source is close 
to the boundary, but they do not easily permit the use of the superposition principle, since the selection 
of complex image locations depends on both the source and the target. We have developed a new, hy- 
brid representation which uses a finite number of real images (dependent only on the source location) 
coupled with a rapidly converging Sommerfeld-like integral. Our method applies in both two and three 
dimensions and its efficiency is illustrated with numerical examples. 



1 Introduction 

A number of problems in acoustics (and electromagnetics) involve the solution of the Helmholtz equation 

(A + k 2 )u(x) = f(x) , (1.1) 

in the half-plane P = {(x, y) E M. 2 ; y > 0} or the half-space S = {(x, y, z) G M. 2 ; z > 0}, subject to suitable 
boundary and radiation conditions. In acoustics, the Helmholtz coefficient k is given by k — — , where u> is 
the governing frequency (assuming time-harmonic motion) and c is the sound speed. In the present paper, we 
assume k is constant throughout the region of interest. For concreteness, suppose that we wish to compute 
the scattered field due to a unit strength point source located at Xq in the presence of a "sound-hard" 
obstacle over an infinite half-space subject to impedance boundary conditions (Fig. [I]). 

We let the total field u tot — u m + u, where u m denotes the (known) incoming field due to the point 
source and u denotes the scattered field. On a sound-hard obstacle with boundary T, the total field must 
satisfy homogeneous Neumann boundary conditions. Since the scattered field involves no sources outside Q, 
it must satisfy the homogeneous Helmholtz equation 

(A + k 2 )u(x) = (1.2) 

for x E P \ in 2D, and for x E S \fl in 3D. On the obstacle boundary T, we have 

du du m 
dn dn 
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Figure 1: Scattering from a "sound-hard" obstacle above an impedance plane. 



Finally, on the interface, we assume the impedance condition takes the form 

du tot 



dn 



iau tot = 0. (1.3) 



There is a substantial literature on impedance problems and we mention only a few relevant papers which 
discuss the computation of the corresponding Green's function [TJ [2j H3 [Til 13 IE] • We assume that Im(a) < 
to ensure that the lower half-space be energy-absorbing J5[ [3] . 

A simple ansatz for solving this scattering problem is to represent the total field as 



i tot (x) = / G imp (x,y)a(y)dy + G imp (x,x ) 



u' 

where Gi mp (x,y) is the Green's function for the half space with homogeneous impedance boundary condi- 
tions. Imposing the Neumann conditions on T yields an integral equation of the form 



~a(x)+ [ d ^(x,y)a(y)d y = - d ^(x,S) 
2 J r on on 



for x G r, where J^- is the outward normal derivative from the obstacle and the integral is interpreted in 
the principal value sense. This integral equation is a Fredholm equation of the second kind, although it has 
spurious resonances for a discrete set of frequencies k. More robust representations are well-known (see [3])- 
In this short note, we restrict our attention to the calculation of the impedance Green's function Gi mp (x, y) 
itself. 

Algorithms for the computation of Gi mp (x,y) date back to the classical work of Sommerfeld, Weyl, and 
Van der Pol pH [13], [12] , who developed both what is now referred to as the Sommerfeld integral and what 
is now referred to as the method of complex images. For a more recent treatment of this problem, see 

[3 m m m mj . 

The main contribution of the present work is the observation that there is an overlooked solution to the 
impedance problem using a hybrid representation involving a finite number of real images which guarantee 
rapid convergence of a residual Sommerfeld- type integral (sections [4] and [5]). Our approach is somewhat 
related to that of Cai and Yu [JJ, which also separates near and far field contributions, but uses an asymptotic 
method for the near field. 

2 Spectral representation of the Green's function 

For k G C with non-negative imaginary part, the solution to 

(A + k 2 )G(x) = —6(xq), (2.1) 



2 



in a homogeneous medium is referred to as the free-space Green's function for the Helmholtz equation, where 
6(xq) represents the usual delta function at Xq. It is well known that the Green's functions are 



G(x,x ) 



-±H {k\x-x \) in 2D, 

ik\x-x \ 

in 3D. 



4tt\x— Xo\ 

A continuous spectral representation of the Green's functions can be obtained by taking the Fourier transform 
of equation (2.1 ). In three dimensions, the Green's function can be written as 

G(x, x ) = - 3 J _ J ^ J ^ X l + X l + Xl -B d\ x d\ y d\ z . (2.2) 

Switching to polar coordinates and evaluating the integral in X z via contour deformation yields the well- 
known plane- wave formula (which itself is sometimes called the Sommerfeld integral) : 

1 r°° e -\/p 2 -k 2 \z~zo\ 
G{X > X0)= 4*L ^k 2 J °MP*P> (2 . 3) 
r = yj{x- x ) 2 + (y- y ) 2 . 

Due to the central role of this formula in scattering theory, there has been much energy devoted to its 
numerical integration. We do not give a comprehensive review of the various schemes available, which are 
largely based on contour deformation, because they are mainly concerned with the square-root singularity in 
the denominator. Our method will address the range of integration required in the Fourier integral parameter 
P- 



3 The impedance problem 



Impedance (or Robin) boundary conditions are used in acoustic scattering to describe surfaces which are nei- 
ther sound- hard (Neumann) or sound-soft (Dirichlet). They are typically derived from the limiting behavior 
of a high-contrast layered medium (where the sound speed is vastly different in the lower half-space). 

If we represent the scattered response to a single point source located at Xq in a form similar to equa- 



tion (2.3|, then it suffices to match Fourier modes between the scattered and incoming solution using equa- 



tion (1.3). We analyze only the 3D case in detail. For this, let the incoming field be given by a point source 
at x = (xo,yo,Zo) located in the upper half-space. The incoming field below x can be written as 



1 f°° e%/ p2 ~ fc2 ( z-2:o ' 
u m (x) = G{x,x ) = — / — J (pr)pdp, 



(3.1) 



where r is given by (2.3). If we assume a similar spectral representation for the scattered field u with 
unknown density a, 



u ( x ) = / / „ ,„ o-(p)J (pr)pdp 1 
Jo VP — * 



(3.2) 



then on the surface z = the impedance boundary condition (1.31 turns into a simple condition on <x(p) for 
each p: 

1 nfl-L*, , s ( 1 e-^ p2 ~ k2zo a(p) 
_ e -Vp -* z « - a(p) - in I 4- — yH ' 

47T 



47T ^/ p 2 _ k 2 yy _ f,2 



0. 



(3.3) 



Solving for the density a, we have 



Hp) 



47T 



-V p 2 -k 2 z 



■\J ' p 2 — k 2 — ia 

Vp 



k 2 



(3.4) 



ia 
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Modulus of the Sommerfeld density 
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Figure 2: For k — 1 and a = 2, the absolute value of the density a is shown for various values of z Q . 



Therefore, the scattered field u can be evaluated via 
u(x) 



1 

47T 



00 e -\/p 2 -k 2 (z+z ) / 



fc 2 



vV - k 2 \ V V 2 - fc 2 + 



J (pr)pdp. 



(3.5) 



Representation (3.5), also called the Sommerfeld integral, can be used to efficiently evaluate the scattered 
field if either z or z are 0(1), since the integrand is exponentially decaying. However, if z + z ~ 0(/i) <C 1, 
then the size of the integration interval must be chose to be 0(ft. _1 ), which can be unreasonably large when, 
for example, a scatterer is near the interface. Figure [2] shows the rate of decay of the density a for varying 
values of zo evaluated at z = 0. 

Since straightforward use of the Sommerfeld integral is clearly not robust for sources and targets near 
the half-space boundary, there have been many attempts at developing more efficient schemes. For example, 
Cai and Yu added a mollifier to the Sommerfeld representation and expanded the remaining high frequency 
components asymptotically [1 . More common is the use of the method of complex images, which we turn 
to next. 



4 The method of images 

Using image sources or charges to impose a given homogeneous boundary condition is a well-known technique 
in classical mathematical physics and potential theory (see [5] for a thorough discussion) . When solving the 
half-space problem with homogeneous Dirichlet boundary conditions, the response to a point source located 
at (x ,yo,zo) is exactly the field generated by a point source of equal and opposite strength located at 
(xn, yo, — z o)- Similarly, for the homogeneous Neumann problem, the response to a point source located 
at (xo,yo,zo) is exactly the field generated by a point source of equal strength located at (xq, yo, — Zq). 
Unfortunately, there is no single image source that can be used in the case of impedance boundary conditions. 
However, it is possible to develop an explicit representation of the impedance Green's function using an 
infinite ray of images, starting at the reflection of the source point (xq, yo, — zq) and continuing vertically 
down (see [TT] for a historical overview). 

Under this assumption, the scattered field u takes the form: 

/>oo 

u(x)= / G(x,x -2zoz-tz)ri(Z)dt, (4.1) 
Jo 
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where z = (0, 0, 1) is the unit normal vector in the z direction. Using the Sommerfeld representation for the 
Green's function G, we then write the scattered field u as 

i r°° f°° e ~\/p 2 -k 2 (z+zo+0 
<x) = — / / === MprHOpdpd^. (4.2) 

In order for this representation for u to satisfy the impedance boundary conditions at z = 0, it is necessary 
that the density r\ satisfy the impedance boundary condition for all p, 



-\/p 2 -k 2 z roo e -VP 2 - fc2 ( z "+'5) / „ f 00 



■q(Od£ + iae~V p2 - k2z ° -ia e^V" 2 ^ 2 ^ ^^)^ = 0, (4.3) 



\J p 2 -k 2 Jq y/ p 2 - k 
which, after some algebra, reduces to a condition on the Laplace transform of n, 

Jo y/p 2 -k 2 +ia , 44 , 

= l-2ia- ' 



VP 5 



fa 



An alternative method for deriving the previous formula is to equate the representation (4.2) with the 
Sommerfeld response formula (3.5 ). In any case, all we now require is the two well-known Laplace transform 
identities: 

-pt* ■ - 1 



£[e"1(») = 



s + (3 (4.5) 
£{S(t-t )}( S ) = e- t ° s , 

where 

POO 

£[/](*) = / e~ st f(t)dt. 
Jo 

Using these formulas, we solve equation ( |4.4| for the density 7j, which is given by 

= 6(0 - 2iae~^. (4.6) 
The full image formula for the scattered field u is then 

poo 

u{x) = G{x, x - 2z z) - 2ia / G{x, x - 2z z - ^e'^d^ . (4.7) 

Jo 

Note that, if a is real, the density r)(£) oscillates throughout its range, so that the decay in the integrand 
comes only from the l/£ decay in the Green's function G. The standard solution to performing numerical 
integration over the infinite interval is to complexify the coordinate £ (see [111 112) ). This is usually done by 
simply replacing £ with — i£ in formula (4.1), yielding 

poo 

u(x) = G(x,x a -2z z)-2ia G(x, x Q - 2z z + i£,z)e~ ai d£ . (4.8) 

Jo 

Instead of an oscillatory density, the complex image contour has enforced exponential decay. The difficulty 
with the formula is that the behavior of the integrand is rather complicated since it involves evaluating the 
Green's function at a complex argument as £ varies. Figure [3] illustrates the norm of the real part of the 
Green's function at two distinct target locations with the same z value. 



5 A hybrid approach 

Rather than using either the Sommerfeld integral or the complex image formula, it is worth reconsidering 
the real image formula, separating it into two parts, a near- field component and a far- field component: 

u(x) = G(x, Xq — 2zqz) 



2ia ( / G{x,x -2z z-Zz)e- ia tdt + J G(x, x a - 2z z - ^)e~ M «d£ 



(5.1) 
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Figure 3: For k — 1 and a = 2, the plot shows the real part of the integrand in the complex image formula, 
and its nearly singular behavior which is dependent on the locations of both the target and the source. 



where C is a parameter of our choosing. From equations (3.5 1 and (5.1 1, it is straightforward to show that 



G(x, x Q - 2z z - fte-^dZ = - / === j— — . M P r)pdp. (5.2) 

4?r Jo VP 2 - k 2 V/ 3 - k 2 + ia 

This is a Sommerfcld-like formula, which for C ~ O(l), decays exponentially fast once p > |fc|, independent 
of z and zo- Therefore, the full impedance Green's function can be written as 

Gim P (x, x ) = G(x, xq) + G(x, x — 2z z) — 2ia / G(x, x — 2z a z — £z)e 

Jo 



lot 
2^ 



oo -y/p^-k^z+zo) -(^/p 2 -k^+ia)C 



(5.3) 



-J {pr)pdp. 



p 1 - k 2 p 2 - k 2 + 

In Fig. |4j we plot the integrand on the right-hand side of (5.2 1 for various values of zo at the interface z = 0. 



6 Conclusions 

We have derived a new formula for the half-space Helmholtz Green's function with impedance boundary 



conditions in three dimensions. The representation (5.3) consists of a short segment of real images and a 



rapidly converging Sommerfeld integral that accounts for the far field. Unlike the method of complex images, 
it permits straightforward use of the principle of superposition. The impedance Green's function is easily 
evaluated to full double precision accuracy, independent of the location of the source and target, using only 
0(k + a) operations and suitable quadratures. 

Although we provided a detailed derivation only for the three dimensional case, the results are nearly 
identical in two dimensions. In two dimensions, repeating the analysis above yields: 



G lmp {x, xq) = G(x, x Q ) + G(x, x a - 2y Q y) - 2ia / G(x, x - 2y y - £y)t 

Jo 



IOL 

27 



oo p -V* 2 -k 2 (y+y ) „-(VA 2 -fc 2 + M )c 



(6.1) 



VX 2 - k 2 VX 2 - k 2 



-dX, 



G 



Difference between Sommerfeld point charge and image chunk 
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Figure 4: For k — 1 and a = 2, the integrand of the far-field Sommerfeld integral (5.2) is plotted for C 
and point charge located at various distances from the half-space boundary. 



where the free space Green's function is now G(x, Xq) = — jHo(k\x — Xo\). 

An obvious extension of our approach is to the evaluation of the layered medium Green's function, 
where each layer has a distinct Helmholtz coefficient and continuity conditions are imposed across the 
layers. Unfortunately, unlike in the impedance case, an explicit image structure is not available. We are 
currently investigating a semi-numerical method that appears promising. Analogous problems appear in 
electromagnetic scattering, and we expect that the hybrid approach discussed here will apply with minor 
modifications. 



Finally, the reader may have noted that formulas (5.1) and (5.3) make use of a discrete image at 
(xo, yo, — zq). As a — > 0, the impedance condition becomes a sound-hard condition and the discrete im- 
age is all that remains, as one would expect. As |a| — > oo, the limit of the impedance condition is a Dirichlet 
condition, while the image formula diverges. It should yield a simple image source at (xo,yo, —Zo) with the 
opposite sign. A complementary image formula can be obtained with the leading discrete image having a 



negative sign by manipulating the Laplace transform. Rather than (4.4), we could write 



OO 







\l ' p 2 — k 2 + ia , 

VS , 6.2 

2Jp 2 - k 2 

= -1+ , V = , 

\J p 2 — k 2 + ia 

and derive a different image formula. In applications a is typically O(l), so we prefer the formula here. 
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